### This script will overlay a density plot of daily max temps (Empirical and AWAP) for a single site

### Load libraries

library('SDMTools')

### Set working directories

out.dir = '/home1/99/jc152199/ChapterOne/Figures/AWAPvsEmpDensity/'

### Read in Raw Data file with preds

raw.data = read.csv('/home1/99/jc152199/brt/FINALOPTIMALMODELS/Model_Data_Plus_Preds.csv',header=T)

### Summarize raw data to monthly values by site
### Start by appending a column for month

raw.data$month = as.numeric(substr(raw.data$date,5,6))
raw.data$year = as.numeric(substr(raw.data$date,1,4))

#### Begin looping by site

for (site in unique(raw.data$site))

	{
	
	### Subset rawdata to a single site
	
	subdata = raw.data[which(raw.data$site==site),]
	
	### Determine density distribution of empirical climate and AWAP
	### Start by removing positions where AWAP is 'NA'
	
	subdata = subdata[which(is.finite(subdata$AWAPmax)==T),]
	
	### Now calculate density
	
	ed = density(subdata$micro_max)
	ad = density(subdata$AWAPmax)
	
	### Start plotting
	
	png(paste(out.dir,site,'.png',sep=''),units='cm',width=18,height=18, res=1000)
	
	### Calculate plotting limits
	
	xlims = range(ed$x,ad$x)
	xlims = round(c(xlims[1]-1,xlims[2]+1),0)
	ylims = c(0,max(c(ed$y,ad$y)))
	ylims = round(c(ylims[1],ylims[2]+.1),1)
	
	### Configure plot space
	
	plot(ed$x,ed$y,type='n',xlim=xlims,ylim=ylims,xlab='Daily Max Temp',ylab='Density',main=site,sub='')
	
	### Overlay density lines
	
	points(ed$x,ed$y,col='blue', type='l')
	points(ad$x,ad$y,col='red', type='l')
	
	### Add a legend
	
	legend('topright',legend=c('Empirical Daily Tmax','AWAP Daily Tmax'), text.col=c('blue','red'), bty='n')
	
	### Shut device
	
	dev.off()
	
	### Show progress
	
	cat('\n',site,'\n')
	
	}
	

